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SUMMARY 


A failure in an aircraft control element must be accommodated in near real time 
by the pilot and/or automatic control system in order to prevent a possible tragic 
accident. A first step in this accommodation is the detection and identification of 
the failure. This report presents the results of an evaluation of the Generalized 
Likelihood Ratio (GLR) technique for the detection and identification of control 
element failures in transport aircraft. 

The GLR technique utilizes the innovations from a Kalman Filter, which estimates 
the aircraft states, as inputs to the GLR algorithm. Under the assumption that the 
failure is a step or impulse vector function in state space, the algorithm computes 
the likelihood of a failure having occurred and an estimate of the failure vector. 

Two such algorithms, one for full GLR and one for constrained GLR, were evaluated 
using a digital computer linear simulation of the longitudinal dynamics of a B-737 
aircraft. The full GLR algorithm assumes that the failure vector can be any vector 
in the state space, while the constrained GLR algorithm assumes that the failure 
vector is from a finite set of vectors in state space. The simulation included 
sensor errors and wind turbulence. 

In the simulation runs all of the hard-over failures were detected by both the 
full GLR algorithm and the constrained GLR algorithm. Soft failures, however, were 
not detected with integration times up to 1.5 seconds. Results of the simulation 
show that while the GLR technique has potential for detecting and identifying air- 
craft control element failures, the effects of wind turbulence on the missed 
detection/false alarm performance and the effects of Kalman Filter model errors are 
significant problems that must be overcome. 


INTRODUCTION 

For certain anticipated failures in transport aircraft operations there are 
established procedures for the pilot to follow. A typical example is the procedure 
for handling an engine outage during takeoff. However, there are nearly an unlimited 
number of additional unanticipated failure modes for which no appropriate emergency 
procedures will be available. These unanticipated failures must be handled by the 
pilot and/or the automatic control system in real time in order to prevent a possible 
tragic accident. 

In the case of a hard-over failure in a control element the pilot may have only 
a matter of seconds to take corrective action before the aircraft reaches an 
irrecoverable condition. In the case of a failure of lesser magnitude the pilot may 
have more time to take corrective action, but the failure and hence the proper cor- 
rective action may be difficult to identify. In either case the pilot may require 
assistance from the aircraft systems to help him take the appropriate corrective 
action in a timely manner. 

A considerable amount of work has been done in the area of failure detection and 
identification (FDI) in dynamic systems, and Willsky has provided a well-known survey 
of many of the available FDI techniques (ref. 1). Chow (ref. 2) and Willsky and Chow 
(ref. 3) have examined the problem of generating residuals from the system 
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measurement data for use in decision making processes to detect and identify 
failures. The detection of failures in sensors has been investigated by several 
authors, including Motyka, et. al. (ref. 4), Deyst, et. al. (ref. 5), Caglayan 
(ref. 6), and Friedland (ref. 7). Beard (ref. 8) and Jones (ref. 9) have developed 
the theory behind the failure detection filter, and Messerole (ref. 10) has applied 
the failure detection filter to the problem of detecting and identifying failures in 
an FI 00 engine. Another FDI technique is the Generalized Likelihood Ratio (GLR). 

The theory of the GLR has been investigated by Willsky and Jones (ref. 11), Chow, 
et. al. (ref. 12), Bueno, et. al. (ref. 13), Bueno (ref. 14), MIT (ref. 15), Liu and 
Jones (ref. 16), and Chang and Dunn (ref. 17). The GLR technique has been developed 
in several versions, including full GLR, constrained GLR, and simplified GLR; these 
versions differ in the a priori assumptions about the type of failure. The technique 
has been exercised in a simplified simulation of the F-8 aircraft dynamics by Bueno 
and others at MIT (refs. 12-15). Tylee (ref. 18) has examined the use of the GLR to 
detect failures in a nuclear reactor. 

This report presents the results of an evaluation of the capabilities of the 
full GLR and constrained GLR techniques for the detection and identification of con- 
trol element failures in a transport aircraft. This evaluation was conducted pri- 
marily by implementing the GLR algorithms in a linear simulation of the longitudinal 
dynamics of a B-737 aircraft and examining the performance of the technique in 
detecting step failures in the elevator, throttle, stabilizer, or spoilers. The text 
includes a brief description of the GLR technique, a presentation and discussion of 
the results, and conclusions. The theory of the full GLR, the theory of the con- 
strained GLR, the aircraft simulation, and the Kalman Filter are discussed in more 
detail in appendices A, B, C, and D, respectively. 
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constrained failure vector 

matrix relating the innovations response to a failure X and defined 

by equations (A34) or (A43) 

filter transfer function 

observation matrix 

altitude, ft 

identity matrix 

index denoting the i-th failure vector from the set of constrained failure 
vectors 

Kalman Filter gain matrix 
sample number 

turbulence scales in the x- and z-directions , respectively, ft 
generalized log likelihood ratio, or log GLR 

•) 
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covariance matrix of the error in the estimate x(k+1 |k) 

A 

covariance matrix of the error in the estimate x(k|k) 
conditional probability density function 

covariance matrix of the plant noise in the Kalman Filter model 

perturbed inertial pitch rate, rad/sec 

pitch rate due to wind, that is, rotation of the atmosphere about the 
y-axis, rad/sec 

covariance matrix of the observation noise v, in the Kalman Filter model 
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covariance matrix of the Kalman Filter plant disturbance ^ 

covariance matrix of the plant noise 

covariance matrix of the wind system plant noise § 
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(•), power spectral density of the x-, angle-of-attack, and pitch rate 
components, respectively, of the wind turbulence 
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time, sec 

time at k-th sample, sec 

trim inertial velocity in x-direction, ft/sec 

control vector 

components of control vector 

steady state wind in the x-direction, ft/sec 

perturbed inertial speed in the x-direction normalized by the trim 
velocity U Q 

normalized speed in the x-direction due to wind, that is, normalized wind 
velocity in the negative x-direction 

gust and steady state components, respectively, of u' 

w 

covariance matrix of the Kalman Filter innovations 
airspeed, ft/sec 

observation, or sensor, noise vector 
wind system state vector 

components of wind system state vector W 
• k 

wind, or plant noise, vector 

steady state wind in the z-direction, ft/sec 

system state vector 

component of system state vector x 

k 

estimate of x at time (sample number) k + 1 given observations through 
time k 

estimate of x at time (sample number) k given observations through 
time k 

Kalman Filter state vector 
observation vector 

component of the observation vector z 

k 

failure magnitude; perturbed inertial angle-of -attack, rad 

part of the angle-of-attack due to winds, rad 

gust and steady state components, respectively, of c^, rad 
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trim value of the angle-of-attack, rad 
vector of innovations from the Kalman Filter 
failure magnitude 
disturbance transition matrix 
Kalman Filter disturbance transition matrix 

disturbance transition matrix for the wind input and defined by 

equation (C55) 

perturbed elevator position, deg 
perturbed stabilizer position, deg 
perturbed spoiler position, deg 
perturbed stabilizer command, deg 
perturbed thrust, klbs 
perturbed throttle command, deg 
unit step function 

Kalman Filter plant disturbance vector 
GLR* threshold 

wind disturbance vector defined by equation (C56) 
failure time, or sample number; perturbed pitch, rad 
trim value of pitch, rad 
likelihood ratio 
failure vector 

components of the failure vector 
wind system plant disturbance vector 

RMS velocities of wind turbulence in the x- and z-directions, respectively, 
ft/sec 

time, sec 

system state transition matrix 
Kalman Filter state transition matrix 





wind system state transition matrix 


control transition matrix 
Kalman Filter control transition matrix 
ft spatial frequency, rad/ft 

OJ , temporal frequency, rad/sec 

Notation: 

~ indicates assumed value or value computed using assumed values for the 

independent parameters 

indicates estimated value 

F superscript F indicates that the variable is used in the Kalman Filter 

formulation 

T superscript T denotes transpose 

E { • } expectation operator 

1 subscript 1 denotes that part of a vector attributed to system dynamics 

with no failure 

2 subscript 2 denotes that part of a vector attributed to a failure 


GLR TECHNIQUE 

As a failure detection and identification technique the Generalized Likelihood 
Ratio method processes the innovations from a Kalman Filter to compute the likelihood 
ratio, that is, the ratio of the conditional probability of the innovations assuming 
a system failure has occurred to the conditional probability of the innovations 
assuming no failure has occurred. This ratio is then compared to a threshold value 
to decide whether or not the system has failed. 

The GLR technique assumes that a Kalman Filter is used to track, or estimate, 
the states of the unfailed system, and that failures in the system are of the type 
that can be described as a vector in state space that occurs as an impulse or a step 
function at time 0. Although the technique is applicable to time-varying systems, a 
time-invariant formulation is used if possible to simplify the computations. Mathe- 
matically, such a system and failure are described by the state equations 
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where 




* 


x^ is the state vector 
is the control vector 

is the vector of disturbances (plant noise) 
z^ is the vector of observations (measurement) 
v^ is the vector of observation noise 
<j> is the state transition matrix 
i|) is the control transition matrix 
T is the disturbance transition matrix 
H is the observation matrix 

A failure in the system is described by the last term in equation (1), where X is a 
constant failure vector which describes the magnitude and direction in state space of 
the failure, 6g ^ +1 is a step (or impulse) function, and 0 is an unknown positive 
integer denoting the time (sample number) of the failure. 

If no failure has occurred, that is k < 0, then from Kalman Filter theory the 
innovations from the filter are a zero mean, white, Gaussian random vector sequence. 
In this case the probability of occurrence of the observations is described simply by 
the multivariate Gaussian density with zero mean. In the case of a failure the 
Kalman Filter ‘innovations are composed of a zero mean, white, Gaussian random vector 
sequence plus a deterministic vector sequence G(k;0)X due to the failure. In this 
case the probability of occurrence of the innovations is described by the multi- 
variate Gaussian density with time-varying mean G(k;0)X. The matrix sequence G(k;0) 
can be computed from the state equations. 

The GLR technique computes the likelihood ratio from the two probabilities just 
discussed and compares the ratio to the threshold to determine if a failure has 
occurred. The GLR algorithms also compute maximum likelihood estimates of the 
failure vector X and the failure time 0. 

The constrained GLR technique differs from the full GLR in the assumptions made 
regarding the failure vector. Whereas the full GLR allows the failure vector X to 
be any vector in the state space, constrained GLR restricts the failure vector to be 
one of a set of vectors f^ with a scalar multiplier, or magnitude, a* This 
technique computes a likelihood ratio for each possible vector f^ for comparison 
with a threshold and computes maximum likelihood estimates of the failure magnitude 
a and failure time 0. 


EVALUATION OF THE GLR TECHNIQUE VIA SIMULATION 
The Aircraft Simulation 

To achieve the objective of obtaining a preliminary evaluation of the capability 
of the failure detection filter to detect and identify failures in an aircraft 
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control system, the filter was implemented in a digital computer linear simulation of 
the longitudinal dynamics of a B-737 aircraft using a small perturbation model. 

These dynamics are described in state space by the equation 


x k+1 = <t* k + + rw k + <*fi + 


(3) 


The aircraft state vector is defined by 
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where 

0 = pitch attitude 

u' = normalized speed in the x-direction 
a = angle-of -attack 
q = pitch rate 
ST = thrust 

SS = stabilizer position 
The control vector is defined by 
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where 


6 e = elevator command (= elevator position) 

6th = throttle command 

6s t = stabilizer command 

6sp = spoiler command (= spoiler position) 

The plant noise vectors and represent the winds and will be discussed 

later. The state transition matrix <|>, the control transition matrix ip, and the 
disturbance transition matrix r were derived from the continuous time system 
matrices A, B, and D, respectively. The system matrices A, B, and D were 
computed from data supplied by the aircraft manufacturer for various aircraft trim 
conditions. All of the results discussed in this report were obtained for the 
aircraft trimmed for final approach on a 3 degree glideslope with the exception of a 
few simulation runs to evaluate the effects of modeling errors. 

To evaluate the performance of the failure detection filter four types of fail- 
ures with three magnitudes for each type were simulated. These four types were step 
failures in the elevator, throttle, stabilizer, and spoiler. The failure vectors 
were derived from the control transition matrix for a step change of unity magnitude, 
for example for an elevator position change of one degree. The three magnitudes were 
chosen to represent a hard-over failure, a soft failure, and a failure of intermedi- 
ate magnitude. The failure types and magnitudes are listed in table I. Of course, 
the failure vectors change as the i|r-matrix changes with the aircraft trim conditions. 


TABLE I.- FAILURE TYPES AND MAGNITUDES 


Failure 

Type 

Magnitude 

Elevator 

Throttle 

Stabilizer 

Spoiler 

10°, 3°, 1° 

40°, 12°, 4°’ 

-6°, 3°, -1° 

8°, 3°, 1° 
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The Wind Simulation 

To provide a realistic wind environment for evaluation of the GLR technique the 
simulation included both turbulence, or gusts, and steady state winds. The gust 
components in normalized x-velocity, angle-of-attack, and pitch rate were modeled 
using the familiar Dryden spectra (ref. 20), which are 
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where 


SI - spatial frequency, ft 

b = wing span 
= 93 ft 

V a = aircraft airspeed, ft/sec 

= 216 ft/sec in the landing configuration simulated 

L U ,L W = turbulence scales in the x- and z-directions , respectively 
= 1750 ft 

o,o = rms velocities of turbulence in the x- and z-directions, respectively, 
ft/sec 

After conversion from the frequency domain to the continuous time domain in 
state space followed by conversion to discrete time, the state equations for the wind 
system become 


W k+1 “ 4 ) W W k + 5k 


(9) 
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where 


W k = 6-dimensional wind state vector with Wy.(5) and W k (6) being the steady 
state components in angle-of -attack and normalized x-velocity, 
respectively 

<j^ = wind state transition matrix 

5^ = zero mean, white, Gaussian random vector sequence 

In order to provide a thorough evaluation of the GLR technique, a wide range of 
wind conditions were included in the simulation runs. RMS turbulence velocities 
corresponding to medium and heavy clear air turbulence (CAT) and to thunderstorm 
conditions were simulated together with steady state winds. Calm (no wind/ 
turbulence) conditions were also simulated. These conditions are summarized in 
table II. 


TABLE II.- WIND CONDITIONS USED IN THE SIMULATION 


Turbulence 

conditions 

RMS 

velocity 

gust 
, ft/sec 

Steady 

wind, 

state 

ft/sec 

°u 

°w 

u g(0) 

CO 

o 

None 

0 

0 

0 

0 

Medium 

2.7 

2.7 

10 

0 

Heavy 

7.0 

7.0 

10 

0 

Thunderstorm 

21 .0 

21 .0 

10,39 

0 


The Measurements 

To enhance the capability to solve the FDI problem it was desired to provide a 
rather complete set of measurements for input to the Kalman Filter and hence to the 
GLR algorithms. On the other hand, for the evaluation results to be credible, the 
measurement set must be technologically feasible, if not typical, for a. modern day 
transport aircraft. The measurements selected for inclusion in the simulation are 
pitch attitude, x- and z-accelerations, pitch rate, airspeed, altitude rate, and 
angle-of-attack. 




The measurements formed a 7-dimensional measurement vector z k defined by 
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The measurements z^d ) - z k (7) include errors such as noise, bias, scale factor, 
and misalignment where appropriate. Numerical values for the errors are listed in 
table C.1 . 


The Filter 

A 10-state Kalman Filter was used to estimate the six aircraft states and the 
four wind gust states and to provide the innovations for input to the GLR algorithms. 
For computational simplicity a constant gain filter was utilized. The measurements 
described in the previous paragraph were not all linear functions of the system (air- 
craft + wind) states as required by the Kalman Filter formulation. To circumvent 
this difficulty a pre-processor was used to transform the mesurement vector into a 
vector of pseudo-measurements, or observations, which could be approximately repre- 
sented as linear functions of the system states. A more detailed description of the 
Kalman Filter and the pre-processor can be found in appendix D. 

The time-invariant full GLR algorithms from equations (A39) through (A52) of 
appendix A were implemented in the simulation for evaluation. At each iteration of 
the simulation the algorithms computed the likelihood ratio ACk;"?) and the estimate 
"XCk; D) of the failure vector for each of the windows of the most recent ? to 30 
innovations. From this set of results the algorithm computed the maximum likelihood 
estimates 'S(k) and X(k) of the failure time and the failure vector, that is, the 
estimates of failure time and failure vector that corresponded to the largest 
likelihood ratio. 

The time-invariant constrained GLR algorithms from appendix B were also imple- 
mented in the simulation for evaluation. At each iteration a series of computations 
similar to those for the full GLR were performed. Estimates a of the failure 
magnitude and 1 of the failure vector index were computed rather than estimates of 
the failure vector. 

The simulation and the GLR algorithm computations were run at a sample rate of 
20 iterations per second. Because the initial interest was in detecting and 
identifying catastrophic type failures, a quick reaction time for the algorithm was 
necessary. Thus, to evaluate algorithm performance only a few seconds of flight were 



required in a simulation run. With these constraints the aircraft could be flown 
open loop with no control system and still not diverge significantly form its nominal 
path before a failure was introduced. Therefore, to simplify the simulation the 
aircraft was flown open loop for all of the results discussed in this report. 

A total of 135 simulation runs were performed to exercise the GLR technique. 

The results of these runs are presented and discussed in the next section. 


RESULTS 
Full GLR 

Filter complexity and Kalman gain .- A series of 33 simulation runs were made to 
evaluate the performance of the GLR technique as a function of Kalman Filter com- 
plexity, or number of filter states. One of the filters was the 10-state filter 
described in detail in appendix D. The second filter was a 12-state filter where the 
two additional states were used to estimate the steady state winds. The third was a 
6-state filter which estimated the aircraft states but no winds. Simulation runs 
were made with each filter for each of the four intermediate level failures with both 
medium and heavy CAT. Runs were also made with no failures introduced. 

The rms estimation errors for the 10- and 12-state filters were similar, while 
the errors for the 6— state filter were somewhat larger. The likelihood ratios for 
the 10- and 12-state filters were very similar, but the GLR's for the the 6-state 
filter were much larger, occasionally by nearly two orders of magnitude, even for the 
no failure case. The larger estimation errors and the very large GLR for the no 
failure case relative to the GLR's for the failure cases resulted in the 6-D filter 
being dropped from consideration even though it was computationally simpler. The 
1 o-D filter was chosen over the 1 2-D because it gave comparable performance and was 
slightly simpler computationally. Furthermore the 12-state system^was unobservable, 
and some problems were experienced with matrix inversion of C(k - 0) because C was 
nearly singular in some cases. The 10— state Kalman Filter and associated GLR 
algorithms then were used to obtain the results to be discussed in the remainder of 
this report. 

The previous runs were made using Kalman Gains computed under the assumption of 
heavy turbulence ( = = 7 ft/sec) for the filter model of plant noise. To 

evaluate the effect of filter gain on the performance of the GLR technique, a 
different matrix of Kalman gains was computed assuming thunderstorm turbulence 
(a _ c = 21 ft/sec) for the plant noise model. Using the new gains five runs were 
then made simulating heavy CAT and intermediate level failures. When compared with 
similar previous runs using the old heavy turbulence gains, the resulting GLR s 
showed only slight differences. Failures were detected in three out of 12 detection 
possibilities with the old gains and in two out of 12 with the new gains. Differ- 
ences were concluded to be insignificant, and the old heavy turbulence filter gains 
were used to obtain the remaining results. 

Thresholds.- Three simulation runs of length 20 seconds each were made under 
conditions of heavy clear air turbulence ( cr = cr^ = 7 ft/sec) with no failures intro- 
duced into the system. Each run was made with a different seed number for the random 
number generator such that each run generated a different sample function, or 
sequence, for the wind turbulence and for the sensor noise. From the total of these 
three runs the largest value computed for the likelihood ratio using a 10— sample 
window of data was selected for the 10— sample threshold value.. Threshold values were 
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similarly obtained from the same runs for windows of 20 and 30 samples. Three 
20 second runs were made rather than a single 60 second run because the aircraft 
system tended to diverge too far from the trim conditions when operating open loop 
for longer times. 

This procedure was repeated under conditions of no turbulence to obtain a set of 
threshold values which could be used in evaluating the effects of wind turbulence on 
the performance of the GLR technique. 

It was realized that this procedure would not establish a statistically reliable 
threshold value suitable for actual aircraf t operation. However, the values thus 1 

obtained were deemed reasonable for purposes of comparing full GLR and constrained 
GLR performance and uncovering potential problems in using the GLR technique for 
aircraft control system FDI . 

Failure detection performance .- Forty-six simulation runs were made to evaluate 
the performance of the full GLR technique under various combinations of failures, 
turbulence standard deviations, and steady state winds. The failures were all single 
point; that is, there were no multiple failures. The likelihood ratios observed for 
each run are shown in tables Ilia, Illb, and IIIc for data windows of 10, 20, and 
30 samples, respectively. Some of these ratios were still increasing when the simu- 
lation was terminated at 5 seconds. 

These ratios were compared with the no turbulence threshold values to determine 
which failures were detected. The results, which are summarized in table IV, show 
that all failures were detected successfully with the exception of the 1 degree 
elevator failure which was detected in only three of the six cases. These results 
are very optimistic in that use of the no turbulence thresholds would produce a 
totally unacceptable false alarm rate. However, the results are useful when compared 
w ^th later results to show the degradation in performance due to wind turbulence. 

The same likelihood ratios were then compared with the thresholds obtained under 
conditions of heavy clear air turbulence, and the resulting detection performance is 
summarized in table V. Under these conditions, that is, when the threshold is set at 
the maximum GLR value obtained in heavy CAT with no failure, all of the hard-over 
failures were detected, regardless of the prevailing turbulence conditions. On the 
other hand, no soft failures were detected, and intermediate failures were detected 
in only eight of 54 opportunities. (There are three detection opportunities in each 
run corresponding to the three data window lengths of 10, 20, and 30 samples.) Eight 
of the runs were repeats of previous runs with new random noise sample functions for 
the turbulence and sensor noise. Examination of table III reveals how the GLR varies 
between the two random sample functions. From table V it can be seen that in three 
cases the failure was detected with one of the sample functions but not the other. 

The 44 instances in which failure detection occurred (not including thunder- 
storms) were examined to determine the time delay between failure occurrence and * 

failure detection, that is, the time that it took after a failure for the GLR to 
cross the threshold. The results, which are summarized in a histogram in figure 1, 
show that half of the detections occurred in less than 1 second. 

Time history plots of the actual system states from the simulation and of the 
estimated states from the Kalman Filter are shown in figures 2(a) and 2(b) for the 
case of no wind turbulence and a -1 degree stabilizer failure. Since there is no 
turbulence and no command input, the aircraft states remain essentially unperturbed 
until the failure at t = 3 seconds. The action of the Kalman Filter in tracking, or 
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estimating, the system statesman also be seen in these curves. Plots of the failure 
vector X and its estimate X produced by the GLR algorithm are shown in 
figure 2(c), and the GLR is plotted in figure 2(d). Introduction of the failure at 
3 seconds does not produce a significant change in the likelihood ratio, and the GLR 
does not reach the threshold and hence does not detect a failure during the 5 second 
run. The GLR is not calculated but is set to zero during the first 2 seconds of each 
run to allow the filter time to settle. 


Similar time history plots for the case of a -6 degree stabilizer failure with 
no turbulence are shown in figure 3. Note the larger perturbations in the system 
states for this case. In figure 3(d) the GLR is seen to cross the thresholds for the 
10-, 20-, and 30-sample windows in the vicinity of 4.2 seconds and thus to declare at 
that point that a failure has occurred. 

Time history plots for the -6 degree stabilizer failure in the presence of heavy 
CAT are shown in figure 4. Notice in this case the significant perturbations of the 
state prior to the failure and the tracking of these perturbations by the Kalman 
Filter. The turbulence and the filter estimate of it can be seen in figure 4(b). 

Once again the plot of the GLR in figure 4(d) shows detection of the failure in the 
vicinity of 4 seconds. Note that the curve of the 30-sample GLR is the smoothest of 
the three as would be expected. 

Failure identification performance.- For the full GLR algorithm failure 
identification consists of estimating the failure vector X. Curves of the failure 
vector components and their estimates for the three stabilizer failure cases in 
figures 2(c), 3(c), and 4(c) show step functions in the stabilizer position X(6) at 
t = 3 seconds, the time of the failures. The other failure vector components,^ X ( 1 ) 
through X(5), are insignificant at the scale of these curves. The estimate X(6) 
is seen to approach the true value; however, in all three cases there are significant 
errors in the .‘estimate X(5). Thus it appears that there are inaccuracies in the 
failure identification results. There are also problems in using the failure identi- 
fication information even if it were accurate; that is to say, it is not clear how an 
estimate X of the failure vector X would be used in restructuring the aircraft 
control system after a failure is declared. 

Model mismatch.- The results discussed to this point were all obtained assuming 
a perfect knowledge of the aircraft system; that is, the aircraft model, or system 
matrices, used in the Kalman Filter computations was identical to the simulated air- 
craft. In actual practice, of course, knowledge of the system will never be perfect. 
To evaluate the effects of inaccuracies in the model on GLR algorithm performance, 
five runs were made in which the Kalman Filter model employed the same landing trim 
conditions used in obtaining the previous results. For the aircraft simulaton, 
however, take-off trim conditions were used. Model differences of this degree would 
be encountered if the same model were employed throughout all phases of the flight. 

The likelihood ratios resulting from these runs are summarized in table VI. For 
each of the four failures the likelihood ratios greatly exceeded the previously 
established thresholds indicating failure detection. The problem, however, is. that 
the likelihood ratios for the no failure run also greatly exceeded the threshold 
resulting in a false alarm. While no attempt was made to determine exactly what type 
or magnitude of model errors could be tolerated without serious degradation in GLR 
performance, this limited amount of data indicates that model inaccuracies are a 
potential problem with the GLR technique. ' 
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Constrained GLR 


Thresholds .- Threshold values for the constrained generalized likelihood ratios 
(CGLR's) were found in a manner similar to those for the full GLR. Three simulation 
runs of length 20 seconds each were made under conditions of heavy CAT with no 
failures introduced into the system. The only significant difference in the pro- 
cedures for the constrained GLR case and the full GLR case was that the constrained 
GLR algorithm for this system produced four likelihood ratios rather than one. 
Therefore, four thresholds were required. The three simulation runs without failures 
were searched to determine the largest value for each of the four CGLR's computed 
with a 10-sample data window. These values then were used as the 10-sampla thres- 
holds for the remaining simulation runs. Thresholds for the 20- and 30-sample data 
window CGLR's were similarly determined from the same runs. 

Failure detection performance .- Forty simulation runs were made with various 
failures, turbulence levels, and steady state winds to evaluate the failure detection 
performance of the constrained GLR technique. The likelihood ratios from these runs 
were compared with the appropriate thresholds to ascertain if a failure had been 
detected. The results, which are summarized in table VII, show that 28 out of 40 
failures were detected. Broken down by length of the data window, the results show 
that 24, 23, and 23 out of 40 failures were detected for the 30- , 20- , and 10-sample 
windows, respectively. In terms of failure magnitudes, all of the 12 hardover 
failures were detected, and none of the eight soft failures were detected. Of the 
20 intermediate failures, 16 were detected by at least one of the data windows. 

The mean detection times (time between failure and first detection) were 1.06, 
0.73, and 0.53 seconds for the 30-, 20-, and 10-sample windows, respectively. A 
histogram of the detection times is shown in figure 5. As might be expected the 
10-sample window tended to produce quicker detection because of the shorter 
integration time. 


Time history plots for the case of a 40 degree throttle failure in the presence 
of heavy CAT are shown in figure 6. The true and estimated aircraft states are 
plotted in figure 6(a), and the increase in thrust due to the throttle failure is 
seen to begin at 3 seconds. In the GLR plot in figure 6(d) it can be seen that the 
failure is detected at approximately 3.5 seconds, or 0.5 seconds after the failure, 
when the GLR corresponding to a throttle failure crosses its threshold. Note that 
none of the other GLR' s cross their respective thresholds until nearly 5 seconds. 

Failure identification performance .- For the Constrained GLR algorithm failure 
identification consisted of estimating i, the failure vector index which defines the 
failure direction, and a, the failure vector magnitude. For the 70 detection 
opportunities in which ^a failure was detected successfully, the failure vector was 
correctly identified (i = i) 52 times for an accuracy of 74 percent. These results 
are summarized in table VIII. Examination of these results shows that the 10-sample 
window did not perform as well as the 20- and 30-sample windows in identifying the 
failure. 

Results of estimating the failure magnitude are summarized in table IX, which 
shows the percent error in the estimate a for the cases where the failure was 
detected and correctly identified. The best accuracy occurred for the elevator 
failure, while the accuracy for the throttle failure was particularly poor. 

For the 40 degree throttle failure discussed previously, examination of the 
estimates of the failure index and magnitude in figure 6(c) reveal that the failure 



is correctly identified as soon as the failure is declared at 3.5 seconds. The 
estimate of the magnitude is approximately 50 degrees at the time of detection, and 
it improves to about 42 degrees after an additional 0.5 seconds. 

Model mismatch .- Eight runs were made using a landing configuration model for 
the filter and a takeoff model for the simulation to assess the effects of model 
errors on CGLR algorithm failure detection performance. Just as was the case with 
the full GLR, the constrained GLR algorithm detected a failure in each of the five 
runs which contained a simulated failure, including one with no turbulence. However, 
even in the three runs with no failure introduced, including one without turbulence, 
one or more of the CGLR's greatly exceeded its threshold resulting in a failure false 
alarm. In fact, in all eight runs a failure was declared prior to 3 seconds. 


CONCLUSIONS 

The application of the Generalized Likelihood Ratio technique to the detection 
and identification of control element failures in transport aircraft has been evalu- 
ated using a linear digital simulation of a B-737 airplane. Hard-over failures in 
the elevator, throttle, stabilizer, and spoiler were successfully detected by both 
the full GLR algorithm and the constrained GLR algorithm. Some intermediate level 
failures were detected, but no soft failures were detected with integration times up 
to 1.5 seconds, that is, with data windows up to 30 samples long. One of the primary 
reasons for the missed detections of the lower level failures was the necessity to 
set the threshold level sufficiently high to avoid false alarms caused by wind turbu- 
lence. Further degradation in performance was caused by mismatch, or inaccuracies, 
in the Kalman Filter model. 

From the results of the simulation runs, it is concluded that: 

1 . The GLR technique has potential application in the detection and identifi- 
cation of aircraft control element failures. 

2. False alarms/missed detections due to wind turbulence and performance degra- 
dation due to system model inaccuracies are significant problems which must be 
overcome before the GLR technique can be used in a practical system. 
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APPENDIX A 


DEVELOPMENT OF FULL GLR EQUATIONS 


The development of the GLR equations in this appendix follows the development in 
reference 11. This appendix is included for the purpose of completeness and to 
include a few details not included in the reference. Furthermore, this development 
will be for step changes, or failures, in the system state, whereas the development 
in reference 11 was for impulse changes. The equations will be developed for a 
general, time-varying, discrete dynamic system and then constrained to a time w 

invariant system. The application to the B-737 aircraft dynamics with its specific 
failures is discussed in the main body of the report and in appendix D. 


The System 

Consider a linear, time-varying, discrete system described by the equations 


x(k+1) = <|>(k+1 ,k)x(k) + i|>(k)u(k) + r(k) w(k) + X6„ 

DfK-rl 


( A1 ) 


z(k+1) = H(k+1 )x(k+1 ) + v(k+1 ) 


( A2) 


where <f>(k+1,k) is the state transition matrix, ijj(k) is the control transition 
matrix, r(k) is the disturbance transition matrix, H(k+1 ) is the observation 
(measurement) matrix, x(k) is an n-dimensional state vector, u(k) is an 
r-dimensional command vector, w(k) is an m-dimensional disturbance (plant noise) 
vector, v(k+1 ) is a p-dimensional vector of sensor noise, and z(k+1) is a 
p-dimensional observation (measurement) vector. The vectors w(k) and v(k) are 
zero mean, white, Gaussian sequences with covariance matrices Q(k) and R(k) , 
respectively, that is 


E |w(k) w T (j)} = Q(k) 6., 

3 K 

E |v(k ) v T (j)} = R(k) 6j k 


(A3) 


( A4 ) 


where is the unit impulse function. 

, A failure in the dynamic system is described by the last term in equation (A1 ) , 
where X is a constant n-dimensional failure vector which describes the magnitude 
and direction (in state space) of the failure, 6 0 k+1 is a unit step function, 
and 0 is an unknown positive integer which denotes the time (sample number) of 
failure. 
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The Kalman Filter 


Now consider a Kalman Filter designed to estimate the state of the system 
described by equations (A1 ) through (A4), assuming no failures. The filter is 
described by 


x(k+l | k) = <|>(k+l ,k)x(k|k) + ip(k)u(k) (a5) 

x(k|k) = x(k|k-1) + K(k) y(Vi) (A6) 

A 1 A I 

where x(k|k) and x(k+1 |k) are the estimates of the state at times k and k+1 , 
respectively, given measurements up to and including time k; <)>(k+1,k) and ijj(k) 
are the state and command transition matrices, respectively; u(k) is the command 
vector; K(k) is the nxp Kalman gain matrix; and y(k) is the p-dimensional vector 
of filter innovations. The innovations are defined by 


y(k) = z(k) - H(k)x (k|k-1) 


(A7) 


According to conventional Kalman Filter theory the variance of the innovations is 
given by 


V(k) = H(k) P (k | k-1 ) H T (k) + R(k) 


(A8) 


and the Kalman gain is given by 


K(k) = P (k | k-1 ) H T (k) V _1 (k) (A9) 

The estimation error covariance is propagated by the equations 

P(k|k) = [I - K(k) H ( k ) ] P ( k | k- 1 ) (A10) 


P(k+1 | k) = <j)(k+1 ,k)P(k k)<j, T (k+1 ,k) + r(k)Q(k) r T (k) 


(All) 


Response to Failure 

Suppose now that a failure occurs at time 9. Since the system is linear, by 
the principle of superposition the system dynamic response can be separated into two 
parts: one part due only to the system dynamics with no failure, denoted by 

subscript 1, and a second part due only to the failure, denoted by subscript 2. 
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(A12) 


x(k) = x 1 (k) + x 2 ( k ) 


z(k) = z 1 (k) + z 2 (k) 


(A13) 


The response of the system to a step failure at time 0 can be found by substituting 
equations (A12) and (A13) into equations (At) and (A2). 


x 2 (6) = X 


x 2 (0+1) = <j>( 0+1 » 0) X + X 


x 2 (0+2) = <(>(0+2,0) X + <)>( 0+1 , 0) X + X 


k 

x (k) = l <)>(j,0)X k > 0 (A14) 

j = 0 


z 2 (k) = H(k)x 2 (k) 


k 

= H(k) l 4>( j , 0) X k > 0 (A1 5 ) 

j=e 


In the same manner the Kalman Filter response can be separated into a part due 
to all effects except the failure and a part due to the effects of the failure. 


A A /S 


x ( k ) 

= x (k) + x 2 ^ 

(A16) 

y(k) 

= y 1 (k) + Y 2 < k > 

(A17) 


The Kalman Filter response to a step failure is more tedious to develop than the 
system response in equations (A14) and (A15). Assume for the moment that the filter 
response can be expressed as 


x 2 (k | k) = F(k; 0) X (A18) 


2.0 


(A19) 


y 2 (k) = G ( k; 0) X 


Expressions for F(k; 9) and G(k; 0) will be developed later. 


The Likelihood Ratio 

The problem now is to determine from an examination of the filter innovations 
whether or not a failure has occurred. This problem can be expressed as an 
hypothesis testing problem, where the hypothesis assumes a failure, that is 


H : k < 0 
o 

H : k > 0 


To decide between the hypotheses, form the likelihood ratio A(k), that is, the 
ratio of the conditional probability of the occurrence of the innovation sequence 
y( 1 ) , y(2), ..., y(k) assuming H 1 to the conditional probability of occurrence of 
the same innovations sequence assuming H Q . The decision will be made by comparing 
the resulting ratio A(k) to an as yet undefined threshold r|. 


H 


1 


> 

A(k) g 


< 

H 


(A20) 


Under the two hypotheses the filter innovations are 

H : y(k) = y (k) (A21 ) 

o 1 

H 1 : y(k) = y (k) + G(k;0)A (A22) 


In order to formulate the likelihood ratio, probability distributions for 0 and 
y must be assumed. An alternate approach is to use the generalized likelihood ratio 
in which the maximum likelihood estimates (MLE's) 0 and \ are used for 0 and 
y, respectively. 

A A 

To find the MLE's 0 and X, we need the conditional probability density 
function for the innovations assuming H, , 0, and X. From the theory of Kalman 

Filters, the innovations sequence y.| ( 1 ) , y-|(2), ..., y^ (k) is a zero mean, white, 
Gaussian random vector sequence. Therefore, the sequence y(1), y(2) , ..., y(k) is 
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a white, Gaussian random vector sequence with means G(1;0)X, G(2;0)X, ...» G(k;0)X, 
respectively, and the corresponding conditional probability density function is 


p[y(1), y(2), y(k)|H 1 , 0, x] 


k 

= n 

j=1 


The MLE X assuming 0=0 is found by finding the maximum of the conditional 
pdf. For convenience use the log pdf, and take the gradient with respect to “X. 


I V< j ) 

( 2 n) 


- 1/2 

P/2~~ 


exp ^ - l [y( j ) - G ( j ; 0) xjV^j) [y( j ) - G(j;0)x] 

L j=i ■ ‘j 


( A23 ) 


V~ (£n{p[y(1), ..., y(k) | h i , 0, X] }) 

1c 

“ - I [~2G T ( j ; 0) V -1 ( j ) y( j ) + 2G T (j;0)v“ 1 (j)G(j;0) X] (A24) 

j“1 


Let the expression in equation (A24) equal zero, and solve for X. 


X(k; 0) = C -1 (k? 0) d(k; 0) 


where 


C ( k ; 0) = l G T ( j ; 0) V _1 (j) G ( j ; 0) 

j = 1 


d(k; 0) = \ G T (j;0) V 1 ( j ) y(j) 

j = 1 


( A25 ) 


(A26) 


(A27) 
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The conditional pdf assuming can now be rewritten as 

yOO |h 1 a] 


k 

n 

j=i 


v( j ) 


1 - 1/2 


(2ir) 


p/2 


exp J - l y T (j ) V _1 ( j ) y( j ) 

l 3=1 


'^r , ^ t 'w aj a/t /v 

+ X C ( k ; 0 ) X - d (k; 0) X - X d(k; 0) 


(A28) 


Similarly, the conditional pdf assuming H Q is 


p[y(1), y(k) | H ] 


k 

n 
j = 1 


|v(j) 


I -1/2 


( 2 tt) 


P/2 


ex p|" I Y T <3) V ’(j) y(j)j 


(A29) 


We are now in a position to write the generalized likelihood ratio A(k), or 
more conveniently the log likelihood ratio £(k; 0, X) as 


£(k; 0, X) 


p[y(1), y(k)|H ,0,X] 

p[y(1 ) , . . . , y(k) | h o ] 


= -^(k;©)! - d T (k; 0 ) x + X^k;©) 


After substituting for X(k; 0) from equation (A25) the log likelihood ratio becomes 


«,(k;0,x) = d T (k;0) C Nk^dCk;©) ( A 30) 

The estimate X(k;0) in equation (A25) and the likelihood ratio £(k;0,X) in 
equation (A30) are ^functions jDf an assumed value 0. We must now find the MLE' 0 by 
maximizing £(k; 0, X) over 0.^ This can be done by computing the estimate X(k; 0) 
and the likelihood ratio £(k;£, X) for all 0 and choosing 0(k) as the value 
of e which maximizes £(k;0,X). 

A 

0(k) = arg max £(k;0,X) (A31 ) 

0 
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and the generalized log likelihood ratio £(k) are then 


The MLE A(k) 


A(k) 

= A(k; 0) 

(A3 2) 

U k) 

- £(k; 0, A) 

A3 3) 


The GLR £(k) in equation (A33) is compared to a threshold, as in equation 
(A20) , to decide if a failure has occurred. If it has, then A(k) and 0(k) are 
the maximum likelihood estimates of the failure vector and the failure time, 
respectively. 

Development of F(k; 0) and G(k; 0) 

We now need to develope algorithms for computing F(k;0) and G(k;0), which 
were previously used but not defined in the equations for the Kalman Filter response. 
Assume that the response at time k-1 in equations (A18) and (Ai9) is correct. Then 
the response at time k is 


x 2 (k|k-1) = <|> < k , k— 1 )F(k-1 ; 0) A 


k 

z (k) - H(k) l <j>( j , 0) A 
j = 0 


k 

Y (k) = H(k) l 4>( j , 0) A - <|>(k, k-1 )F(k-1 ; 0) A 

j=e 


= G(k; 0) A 


k 

. * .G(k; 0) = H(k) J <j>(j,0) - <J>(k,k-1 ) F(k-1 ; 0) k > 0 (A34) 

j=0 


x 2 (k|k) = 4»(k,k-1 ) F(k-1 ; 0) A + K(k)G(k; 0) A 


• * .F(k; 0) = (Ji ( k , k- 1 ) F ( k- 1 ; 0 ) + K(k)G(k; 0) k > 0 (A35) 


Equations (A34) and (A35) are recursive relationships for G(k;0) and F(k;0). 
A starting condition for computation can be found by considering the filter response 
at time k = 0. A failure at time 0 has no effect on the predicted response 
x( 0| 0-1) . Therefore 
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A 

x (e|e-n = o 
z 2 ( e) = h( 0) x 
y 2 (0) = H( 0) x 
.* .G ( 0; 0) = H( 0) 

X 2 ( 0 1 0) = K( 0)H( 0)X 
.* .F( 0; 0) = K( 0 )H( 0 ) 


A36) 


(A37) 


For k < 0, the failure has no effect, and 


G(k; 0) = F(k; 0) = 0 


k < 0 (A38) 


Time Invariant Systems 

The algorithms in equations (A25) through (A27) and (A30) through (A38) consti- 
tute quite a computational workload for real time operation. This workload can be 
reduced significantly for the case of a time invariant system and a constant gain 
Kalman Filter.* From equations (A34) and (A35) the filter response for this case is 


G(k; 0) = G(k-0) 

k < 0 

k = 0 (A39) 

k-0 

l <|>(j) - <(>(1 ) F(k-1 - 0) k > 0 

j=0 

F ( k; 0) = F ( k- 0) 

r 

0 k < 0 

_ J k ™ 0 ( A40 ) 

4>( 1 )F(k-1 -0) + KG (k-0) k > 0 

V 


r 

0 

= 1 ” 
H 
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The matrices C(k;l$) and d(k; 0) now become 
C(k; 0) = C(k- 0) 


G T (j-0) 


G ( j-0) 


The sequences of matrices G(k-K), F(k-0), and C(k-0) can all be calculated off- 

line and stored. The matrix sequence d(k? 0) must be calculated in real time since 
it depends on the sequence of innovations. 

Although the computational burden is lessened for a time invariant system, it is 
still formidable. The biggest reason for this is that the formulation presented 
above calls for the GLR and the estimates to be computed at each time k for all 
possible failure times 0 from zero to the present (time k). Thus, the number of 
computations increases without limit as time passes. To limit the computational load 
the algorithms can be defined to use only a window of data at any time; that is, only 
the N most recent innovations ) t y(k-N+2) , . .., are used in the 

calculations. The algorithms for the off-line quantities now become 


k 

= l 

(S# 

j-9 

d(k;0) = \ 

IV 

j=e 


— 1 /v 

V G( j-0) 


V -1 y(j) 


(A41 ) 


( A42 ) 


G(k-0) = G(n) 

H 

— < n 

H l <f>( j ) - fd ) F(n-1 ) 
^ j=0 

F(k-'S) = F(n) 
r 

KH 

< 

if>( 1 )F(n-1 ) + KG(n) 

V 

C(k-'g) = C(n) 

n 1 

= l G (i) V~ G(i) 
i=o 


0 < n ^ N-1 
n = 0 

n > 0 ( A43 ) 

0 < n < N-1 

n = 0 (A44) 

n > 0 

0 < n < N-1 

( A45 ) 
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The computations that are done each sample time k become 

k 

d(k;6) » l G (j-0) V ^(j) k-N+1 < 0 < k (A46) 

j=e 

A(k;0) = C _1 (k-0)d(k; 0) k-N+1 < 0 < k (A47) 

£(k;0) = Y^k; 0)C(k-0) - d T (k;0) A(k; 0) 

- X* 1 Oc; 0)d(k; 0) k-N+1 < 0 < k (A48) 

A 

0{k) = arg max £(k;0) k-N+1 < 0 < k (A49) 

0 

A(k) = X(k;0) A50) 

£(k) = Jt(k; 0) (A51 ) 

A failure is declared if £(k) exceeds a threshold, that is 

H 

1 

> 

n (A52 ) 

< 

H 

o 
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APPENDIX B 


DEVELOPMENT OP ALGORITHMS FOR CONSTRAINED GLR 


The development of the algorithms for constrained GLR is very similar to the 
development for the full GLR in appendix A. In this appendix we will concentrate on 
the differences between the two. 

In the theory of the full GLR the failure vector X can be any vector in the 
state space. However, in constrained GLR theory it is hypothesized that only a set 
of L failure vectors f^ are possible. Thus the state equations become 

x(k+1) = i()(k+1 ,k)x(k) + ip( k)u(k) + r(k)w(k) 

+ ocf i <S 0fk+1 i -'1,2, ..., L <B1) 


z(k+1) = H(k+1)x(k+1) + v(k+1 ) (b2) 


where the f^ are n-vectors representing L failure directions and a is a scalar 
representing the failure magnitude. The other quantities are the same as for the 
full GLR. 


. Filter Response 

A Kalman Filter, designed for the no-failure case, is used to estimate the 
system state. In the case of a step failure with failure vector f^ at time 0, the 
system and filter responses can again be separated into two parts: a response due to 
all effects other than the failure and a response due to the failure. The system 
response to the failure is 


k 

x 2 (k) = l <|>( j , 8)f . a 

j=e 1 


i — 1, ..., L 
k > 0 


k 

zfk) = H(k) l <},(j, 0) f. a 

2 j = 0 


i — 1, ..., L 

k > 0 


The response of the Kalman Filter to the failure is 


x 2 (k|k) = F(k; 0) f<a 


i — 1, »««, L 


(B3) 


(B4 ) 


(B5) 


y (k) = G(k; 0) f . a 

2 l 


i — 1, ..., L 


(B6 ) 
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where the matrices F(k; 0) and G(k; 0) are computed as in equations (A35) and 
(A34), respectively* of appendix A. 

The conditional pdf for the innovations sequence y(2), y(k) 

assuming (assuming a failure has occurred) is 


p[y(1),..., YOO |h i ,0,i,a] 


k 

n 

j-1 


iv(i>r ,/2 


( 2n) 


p/2 


exp < - £ [y(j) ~ G(j;0) fr a] 

L j=1 


. V \j) [ Y ( j ) - G(j; 0) fr ajj 


The MLE a of the failure magnitude assumming , 0, and i is obtained 

by finding the maximum of the conditional pdf. 


(B7) 


— - [y( 1)/ •••/ y(k)| H ,0,i,a] = 0 
3ct 


Solving for a gives 


fr T d (k; 0) 

a(k; 6,i) = — 

fr C(k;0) fr 


(B8 ) 


where C(k;0) and d(k; 6) are defined in equations (A26) and (A27), respectively, 
of appendix A. 


The Likelihood Ratio 

The generalized log likelihood ratio can be written as 


£(k;0,i,a) = d T (k;0) 


fr 

i 


/•sj /N/ T 

a + a f~ d(k; 0) 
1 


rss T* rsj 

a f~ C ( k ; 0) 

l 


f~ 

l 


a 


T rv T <v 

d (k; 0) fr fr d (k; 0) 

1 1 

fr T C (k;0) fr 
l l 


(B9) 
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Similarly to the procedure for full GLR, the estimate a(k;0,I) and the 
likelihood ratio &(k? 0,T, a) are calculated for each 1 and 0, and the MLE's of 
0, i, and a are those which maximize £(k; 0, i, a) . 


0,i = args max Jl(k;0,i,a) 
6,i 


a(k) = a(k; 0,i, ) 


£(k) = £(k;0,i,a) 


(BIO) 


(B11) 


(B12) 


Just as in the case of full GLR, the computations for the constrained GLR 
simplify when the system is time invariant. Similarly, the computational workload is 
further reduced by using only a window of thfe N most recent measurements. The 
development of the algorithms for these more restrictive cases proceeds very 
similarly to the development in appendix A and will not be repeated here. 
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APPENDIX C 


SIMULATION DESCRIPTION 


The Aircraft 

The system used to evaluate the GLR techniques was a discrete, linear, small 
perturbation simulation of the longitudinal channel of a B-737 aircraft similar to 
the simulation described by Halyo in reference 19. Much of the development in this 
appendix follows that by Halyo. The discrete system was derived from a continuous 
time system described by the following state equations: 


x(t)=Ax(t)+Bu(t)+Dw(t) (Cl) 

where the state vector x(t) is a 6-component vector defined as' 


0 

u' 


x(t) = 


a 

q 

<5T 

fiS 


(C2 ) 


and where the perturbed states are 
0 = pitch 

u' = normalized x-velocity 
a = angle-of-attack 
q = pitch rate 
fir = thrust 

fiS = stabilizer deflection 

The thrust and stabilizer states were included to account for the engine spool 
up/spool down time and for the time constant in the stabilizer actuator. The command 
vector u ( t ) is defined by 


u ( t ) 


<5e 

6th 

fist 

fisp 


(C3 ) 
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where 


6e = elevator command (= elevator position) 
6th = throttle command 
6st = stabilizer command 

6sp = spoiler command (= spoiler position) 
The wind vector is defined by 


w( t) 


u' 

w 

“w 



(C4 ) 


where the wind components are 

u^ = normalized wind velocity in the negative x-direction 

c ^ = part of the angle-of-attack due to winds 

q w = rotation of the atmosphere about the y-axis 

Unless otherwise specified, the quantities are defined and the equations are written 
in the aircraft's stability axes. 

Values for the system matrices were obtained using a computer program (TCVOPL) 
which computes the aerodynamic coefficients for aircraft trim conditions specified by 
the program user. Exceptions to this procedure are the coefficients for the thrust 
and stabilizer states. The engine thrust is modeled as 


6T = 


A 55 6T + B 52 6th 


= -0.5 5T + 0.298 6th 


(C5) 


The value 0.5 was approximated from engine data for the B-737. 

The response of the stabilizer on the B-737 is very slow. In order to have 
another control surface for restructurable controls, the time constant of the stabi- 
lizer was artificially shortened to make the surface useful. The stabilizer dynamics 
were assumed to be 


6 s “ ^6 53 + B 63 fist 


(C6) 


= -0.667 6S + 0.667 6st 
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The Winds 


The wind effects on the aircraft are simulated by adding them to the other aero- 
dynamic forces acting on the aircraft. As noted in equation (C4) the wind is com- 
posed of x-velocity, angle-of-attack, and pitch rate components, and they include 
both steady state winds and wind gusts, or turbulence. 

The gust components are modeled using the familiar Dryden spectra (ref. 20), 
which are defined by the following: 


S u (fl) = 


2a 2 L 
u u 


i + (l nr 
u 


S ja) = 


ct 2 l i + 3(l n) 2 
w w w 


2 ^2 


V* [l + (L fi) 1 
a L w J 


f^V 2 

s 0 (n) = - — r • s (n) 

3 „ 2 a 


1 + 


4bn 


s q« (B) 


(C7 ) 


(C8 ) 


(C9) 


(CIO) 


2 2 

In equations (C7) through (CIO) a and a are the variances of the gust veloci- 
. . u w 

ties m the x- and z-axes, respectively, L u and L w are the turbulence scales in 

these axes, b is the aircraft wingspan, V a is the aircraft airspeed, and ft is 
the spatial frequency of the turbulence, which is related to the temporal frequency 
to by 


to = fiV a 


(C11 ) 


In order to use these spectra in the simulation, which generates random gusts as a 
function of time, the spectra must be converted from functions of spatial frequency 
to functions of temporal frequency using the relationship 


= ~~ sU/v a ) 

a 


S( o>) 


(Cl 2 ) 



Angle-of -attack component .- Upon conversion to a function of 
density of cig becomes 


10 


the spectral 


S a ( <o) 


a 2 L 1 + 3(L /V ) 2 w 2 
w w w a 

[l + (L /V ) 2 a> 2 ] 2 
a 1 w a 1 


(Cl 3) 


which can be factored into 


S Q ( a>) 


a 2 L 
w w 


1 + jV3 (L /V )io 1 - j/3 (L /V )<o 
w a w a 

• • " — - — - 

[! + j(L w /V a )o,] 2 [1 - j(L w A a )u,] 2 


(Cl 4) 


Define a transfer function G (s) using the realizable part of the spectrum in 
(Cl 4 ) . 


G a ( S ) 


1 + /3 (L /V )s 
w a 

[i + ayv a ) S ] 2 


1 + /3 (L A )s 
w a 

1 + 2 (L A )s + (L A ) 2 s 2 

w a w a 


(Cl 5 ) 


A filter with this transfer function driven by white noise with a density 
will produce random gusts with the spectral density specified by equation 
with variance 


2 

of a L /V 
(Cl 3) W and 


3 

a 


Let us now turn our attention to obtaining a set of state equations which de- 
scribe the filter specified by equations (Cl 5). First convert the transfer function 
to an equivalent scalar differential equation. 


'L \ 2 2L 

w \ w * 

v~ °g + v~ a g + °g 
a / a 


/3 L 
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2V 

•• 3 • 

+ r~ °g + 




(Cl 6) 
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Now let y = - 0 g and ^ = -£ in equation (C16), and let 
w i = y = -Og 

W 1 = w 2 + Cl 

W 2 = -2(V a /L w )W 2 - (V a /L w ) 2 Wl + c 2 ^ 

From equations (C17) find expressions for y and y. 

y » w 2 + c i ^ 

y - -2(V a /L w )W 2 - (V a /L w ) W 1 + c 2 £.j + ^ 

Substitute these in equation (C16), and solve for and c 2 

Cl ^ * (V a /L w ) 

c 2 = (1 - Si2)(V a /L w ) 2 


(C17) 


(Cl 8 ) 


(Cl 9 ) 


Pitch component .- Upon conversion to a function of oj the power spectral 
density of the pitch gusts becomes 


S q U) = 


0) 


1 + (4b/W )V “ 

a 


S (o>) 


(C20) 


1 + j(4b/nv)a> G a (ja)) * 1 - j(4b/irV ) u G a ( “ ja)) 


-1 QJ 
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The wind gusts in pitch can be simulated by passing angle-of-attack gusts, that is 
the output of the filter G a (s), through a filter with a transfer function Gg(s) 
defined by 


G q (s) 1 + ( 4b/nV ) s 

3 


(C21 ) 


To obtain a state variable representation of the filter, first find the scalar 
differential equation equivalent of the transfer function Gg(s). 


4b 

uV 


a 


q g + 


% = 



(C22) 


Let 


q 


g 


irV 
a 

4b~ 


W. 


uV 

C 

4b" 


W. 


(C23 ) 


Then 


. * • * V a * V a 

- w l + W 3 - 4b~ W l + W" W 3 


= -W, 


(C24 ) 


Therefore 


ttV ttV 

• a a 

W t - ^TT~ W_ 

3 4b 1 4b 3 


(C25 ) 


X-axis component .- Upon conversion to a function of to the power spectral for 
the gusts along tihe x-axis becomes 


S u (co) 


or 2 
2L a 
u u 


1 


1 + (L /V )V 

u' a 


(C26 ) 


2L a 2 
u u 


1 + j(L /V ) to 1 
J u a 


J‘W" 
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2 3 

The x-axis gusts can be obtained by passing white noise with density 2L a /V 
through a filter with transfer function G u (s) defined by 


G u (s) 


1 

1 + (L /V )s 
u' a 


(C27) 


The corresponding scalar differential equation is 


V 


i' 

g 


(C28) 


or in terms of state variables 




(C29 ) 


Steady state winds .- The steady state winds have an x-axis component and an 
angle-of-attack component each of which is simulated as the output of a first order 
differential equation driven by white noise with a very small variance. This allows 
the winds to vary slightly during a run. The desired steady state wind velocities 
are used as the initial conditions. Mathematically the winds are described as 
follows: 


u' (t) = F ; u' (0) = 
s 3 s 


a (t) = F ; a (0) = 
s 4 s 


u (0)/U 
s o 


w (0)/U 
s o 


(C30 ) 


Wind state equations .- Equations (Cl 7), (C25), (C29), and (C30) can be combined 
to describe the total winds as follows: 


W = V* + 


(C31 ) 
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The covariance matrix of the Gaussian vector £ is 


F = eU5 T } 


F 



(C35 ) 


In the foregoing formulation the wind gusts W 1 through W 4 were written in 
aircraft body axes, and the steady state winds VJg and Wg were written in earth 
axes. They must be transformed into stability axes. Furthermore, the gust and 
steady state components must be added, and and Wg must be combined according 

to equation (C2) in order to obtain proper wind forces acting on the aircraft. This 
can all be accomplished with the transformation as follows: 


w = C^W 


(C36) 


where 
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(C37 ) 
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Discrete State Equations 

The state equations for the aircraft and winds have been written as functions of 
continuous time t. These equations need to be coverted to functions of discrete 
time t k for use in a digital computer simulation. Such a discretization has bden 
discussed by Halyo (ref. 19) . 

Consider the aircraft continuous time state equation (Cl ) . If the state is 
known at time t k , the state ate time t k+/ j = t k +■ T can, be found by integration over 
the interval t k to t k+1 . The result is shown in many texts (e.g. , ref. 21) to be 


x(t k+i ) = ^k+i ' t k )x(t k ) + / 4 ’ <t k+r T)Bu(T) dT 


t; 


- I 

t. 


'k+1 


k+1 


, t)Dw( t) dx 


where 


<j)(t, x) = e A(t " T) 



(C39 ) 


To simplify the notation, let x(t k ) be denoted by x k ,, and assume that the command 
u(t) is constant over the integral t k < t < t k+1 . Then 


k+1 

x k+1 = <M k + 1 ' k)x k + / ^ (t k+1 ' t) dT Bu k 


+ 1 

+ / (}i(t k+i , x )Dw( x) d- 

t k 


{ C40 ) 


The integral 


/ k+1 * (t k + l' T) dT 


can be expressed as 


^k+l 

J e 


A(t k + r T) 


dx B = e Atk+1 / tk+1 e' AT dx B 


(C41 ) 
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If A has an inverse A 1 , then 


* . J tk+1 e Mt ^r X) dT B = / tk+1 e- a V' <3 ( At)b 


At. 


k+1 


-At 


-e 


k+1 


' N ,- 1 

+ e I A 
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(C42) 


= [ 4>(k + 1 ,k) - I]A -1 B 


The aircraft state equation now becomes 


k+1 


k+1 

<()X k + + / (|)(t k+1 , t)Dw( t) dx 

fc k 


(C43) 


To evaluate the integral involving the wind w(x) r we first convert the 
continuous time wind state equation (C31 ) to discrete time. 


W = *A W W + B w 5 


(C31 ) 


Upon integration this becomes 


w k+1 = <t>W W k + ^k 


(C44) 


where ^ is the state transition matrix for the wind system. 




WrV 

— e 


(C45 ) 
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and 


Ek =■ T :Jt * 

% 

(C46 ) 

. /■ ■V T ’-->\E'- k * •> 

;0 


Halyo (ref. 19) has shown that is a Gaussian white noise vector sequence with 

covariance Rg. 

R 5 (k,j) - Ef;^) 



f ‘ K (T '’ ! »wE !t k 


+ t) dx 




T T 

J J 4) w (T,T)B w E{5(t k + T)C (t. + s)} B w «j, w (T,s) dx ds 


Since £(t) is a white noise process, and since the intervals ( , t^ + T) and 
(tj'tj + T ) do not overlap for j*k. 


E { £ ( t^ + x)g T (tj + s) } 



j * k 
j = k 


(C48 ) 


Therefore, 


R £(k, j ) = R^fijk (C49) 

T 

H = / y T ' s)B W PB wV T " S) ds ' (C50) 
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The matrix C^ t which transforms the continuous time wind variables W into the 
variables w for use as the plant noise in the aircraft state equations, is still 
valid for the discrete case. Thus 


w k = Cw W k 


(C51) 


We now return to the evaluation of the integral in equation (C43) involving the 
wind w( x) . A change of variables produces 




<{>(t k+1 , x)Dw( t) dx = 


T 

j <J>(T,s)Dw( t^ + s) ds 
0 K 


(C52) 


After substitution of equation (C36), integration of equation (C44) from 0 to T, and 
substitution of equation (C46), the integral becomes 


T 

/ <j>( T, s )Dw( t, +s) ds 
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T 


I 

0 


ds w k 


T S 

+ / <j>(T,s)DC / <fc (s,x)B £(t + x) dx ds 

0 0 

(C53) 


The aircraft state equation can now be written as 


x k+1 = ^k + ^ u k + r W w k + r 'k 


(C54 ) 


where 


r w = / <J)(T,s)DC w (( lw (s,0) ds 

0 


(C55) 


and 


T 


hjc = / <}>(T f s )DC / 4^(s,x)B „?(t + x) dx ds 

0 0 


(C56 ) 
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The sequence is a Gaussian white noise vector sequence with covariance R 


V k ' j > = 


T T 

= / / <(,(T f s)DC w <^(s,T) ds B w 5(t k + T ) dx 

Ox 



T r T T T T 

B w / <f> w (x,y)C w D <|> (T,x) dx dy dx 


As we noted previously in equation (C48), 
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(C58) 


(C59) 


• / ^(x,x)cy> T (T,x) dx 

x 


dx 


(C60 ) 


This completes the discretization of the state equations necessary to simulate 
the aircraft and winds. 
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Measurements 


Thus far we have described the simulation of the aircraft dynamics and of the 
winds acting on the aircraft. For the purposes of the work discussed in this report 
the simulation needs one additional element - simulation of the measurements, or air- 
craft sensor outputs. 

The measurements of interest are pitch attitude from a pitch gyro, x- and 
z-accelerations from body mounted accelerometers, pitch rate from a pitch rate gyro, 
airspeed and altitude rate from the air data computer, and angle-of-attack from an 
angle-of-attack vane. A reasonably accurate simulation requires computation of the 
true value of the quantity being measured and then the addition of appropriate 
errors. For each measurement these errors include some, but not necessarily all, of 
the following: bias, white noise, scale factor, and alignment. The true values and 
errors are combined at each sample time t k to form a measurement vector z k . 


z k < 1 > 

z k (2) 



(C61 ) 


Pitch.- The pitch measurement includes additive noise and bias errors. 


z k (1 ) = 0 + v k (1 ) + b(1 ) 


(C62 ) 


where 


0 = pitch = x k (1) 
b( 1 ) = 0.23° 

v k (i) = zero mean white Gaussian noise 
e { v £( 1 * 1 " <°‘ 23 °> 2 

Accelerations .- The x- and z-axis acceleration measurements include noise, bias, 
scale factor, and misalignment errors. The first step in obtaining a simulated mea- 
surement is to compute the true value of the acceleration in stability coordinates, 
A x and A z , from the equations of motion. 
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(C63) 


Ax g = U 0 u' + qu Q tan o + g sin( 0 O + 0 - a Q ) 


A z = *oP tan a + d + u') a/cos 2 a] - qU Q - g cos(0 Q +0- a Q ) 


where 


g = acceleration due to gravity 
U Q = trim inertial velocity along x 
0 Q = trim pitch attitude 
Oq = trim angle-of -attack 
0 = perturbed pitch = x k (1) 

u' = perturbed inertial velocity along x = x k (2) 
a = perturbed inertial angle of attack = x k (3) 
q = perturbed pitch rate = x k (4) 

At any time t k values of x k are known from the aircraft state equation. The 
quantities u' and a are calculated using the coefficients from the continuous 
time state equation. 


u' = £ ^[Ax k + Bu k + Dw k ] 
a = ^[Ax k + + DwJ 


> (C64) 


where 


£^=[0 1 0 0 0 0 ] 
£^=[0 0 1 0 0 0 ] 
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These stability axis accelerations are transformed into body axis accelerations, and 
the effects of accelerometer misalignment are added as follows: 


A v = A costa- + a) - A sinta.. + a) 
x b ^s ° z s ° 


(C65) 


A z ^ = A x sin(c^ + a) + A z costc^ + a) 


J 


where 

a = accelerometer alignment error = 0.2° 

Noise, bias, and scale factor errors are added to obtain the acceleration 
measurements. 


z k (2) 
z k ( 3) 



+ Vv. ( 2 ) + b( 2 ) 


I 


+ v k ( 3 ) + b ( 3 ) 


'] 


(1 + s) 
(1 + s) 


j (C66) 


where 

s = scale ‘factor error = .0025 
b ( 2 ) = b(3) = bias = 0.32 ft/sec 2 
v k (2) , v k (3) = zero mean Gaussian noise 
e{v 2 (2)} = e{v 2 (3)} = (0.32 ft/sec 2 ) 2 

Pitch rate.- The pitch rate measurement includes only a noise error. 


z k (4) = q + v k (4) 


(C67 ) 


where 

q = pitch rate = x k (4) 

v k (4) = zero mean Gaussian noise 

e{v 2 (4)} = (.02 deg/sec) 2 
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Air speed .- The air speed measurement includes a multiplicative noise error and 
a bias. 


.(5) « V [1 + v (5) ] + b(5) 


(C68 ) 


where 


V a = U Q (1 + u' + u^)/cos(a + o^) 

u' = normalized inertial velocity along x = x k (2) 
u^ = normalized wind velocity along x = w k (1) 
a = angle-of-attack due to inertial velocity along z = x k (3) 
cfy = angle-of-attack due to wind velocity along z = w k (2) 
b(5) = 3 kts 

v k (5) = zero mean Gaussian noise 
e{v 2 (5) 1 = (0.02) 2 

Altitutde rate .- The altitude rate measurement includes only a noise error. 

z k C6) = h + v k (6) (C69) 

where 


h = U Q (1 + u' ) [sin( 0 Q + 0 - a Q ) - tan a cos(0 Q + 0 - a Q ) ] 
v k(6) = zero mean Gaussian noise 
e{v 2 (6) } = (5 ft/sec) 2 

Angle-of-attack .- The angle-of-attack measurement includes bias and additive 
noise errors. 


z k (7) = a+ + v k (7) + b(7) 


(C70 ) 


where a and were previously defined, and v k (7) is zero mean Gaussian noise. 

The bias b(7) and noise error variance were estimated to be 0.25° and (0.4° ) , 
respectively. 
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The measurement errors are summarized in table Cl 


TABLE Cl MEASUREMENT ERRORS 




Error type 


Measurement 

Noise 

standard 

deviation 

Bias 

Scale 

factor 

Alignment 

Pitch 

0.23° 

0.23° 

- 

- 

Acceleration 

0.32 ft 

2 

sec 

0.32 ft 

2 

sec 

0.0025 

0.2° 

Pitch rate 

0.02 deg 
sec 

- 

- 

- 

Air speed 

0.02* 

3.0 kts 

- 

- 

Altitude rate 

5 

sec 

- 

- 

- 

Angle-of-attack 

o 

• 

o 

0.25° 

- 

- 


♦multiplicative 


. Implementation 

Implementation of these equations into the simulation requires evaluation of the 
integrals in the expressions for r w , R^, and R^ in equations (C55), (C50), and 
(C60), respectively. The integrals in equations (C55) and (C50) and the inner inte- 
grals in equation (C60) were evaluated using a Langley software library subroutine 
GLEGEN, which performs numerical integration using a Gauss-Legendre formula. The 
outer integral in equation (C60) was evaluated using the library subroutine SIMP, 
which performs numerical integration using Simpson's formula. In all cases, the 
aircraft transition matrix was evaluated using the library subroutine CONEXP, which 
computes the matrix exponential. 

Random sequences .- The random sequences 5^ and have correlation matrices 

Rj- and R^ defined by equations (C50) and (C60), respectively. In general, these 
matrices are not diagonal, and thus the components 5^(1), £^(2), ..., £^(6) of the 
vector are not independent as they were with the vector £(t) in the continuous 

time case. This is also true of the components 11^(2), ..., ^(6) of n^. 

At time tj^ the simulation generates a vector of six random numbers with 

covariance matrix R^ using the following technique. Let x be a vector of zero 
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mean, independent, Gaussian random variables with unity variance. Then the correla- 
tion matrix R x is 


T 

R x = E{xx } 
= I 


(C71 ) 


Let y be the vector of desired zero mean random variables y 1 , y 2 , y N with 

the covariance matrix Ry, Now let the desired random vector y be given by 


y = Gx (C72) 

where the transformation G is defined to produce the desired covariance of y, that 
is 


T T T 

E{yy } = E {Gxx G } 


m m 

= GE{XX }G 


(C73) 


GG 


T 


Therefore 


GG T = Ry (C74 ) 

If G is assumed to be triangular, one solution (of many) to equation (C74) is given 
by the following: 

Y 1 = g 11 x 1 
y 2 = g 21 x 1 + g 22 x 2 

• (C75 ) 


y N - g N1 x 1 + g N2 x 2 + ••• + 9 nn x N 
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Then 


E{y^} - R,, - g*, 

E(y 2 y, I - » 21 - 9„g 21 

E {y 2 } = R 22 - g 21 + g 22 (C76) 


This system of equations can be solved for the elements 


’ij 


to give 



g j1 “ nj 1 / g 1 1 


2 < j < N 


j i 



i-1 

I g. g. 

im im 

m=1 



< i < j 


(C77 ) 



j-1 


, 1/2 


V 2 

l ^ 


m=1 


3 m 


2 < j < N 


The same technique is used to generate the vector of random numbers with 

covariance R^. 

Failures Failures in the control elements were simulated as steps or ramps in 
the state variables according to 


x k+1 = + ^k + r W w k + % + a k f i 6 k+1,0 


(C78 ) 


51 



where 


f i = 


h 


= i-th column of the control transition matrix, 
ahd is chosen to be a step or ramp. 

Because the acceleration measurements z^d) and z^(3) contain terms which 
ihclude the continuous time control input matrix B as in equations (C64) through 
(C66), any failure affecting these terms must be accounted for in the simulated mea- w 

surements. For the four failures simulated, only the elevator arid spoiler failures 
affected the B-matrix and thus introduced failure effects directly into the measure- 
ment equation via the expressions for u' and a. For a stuck actuator these 
effects were simulated by adding terms to equation (C64) as follows: 


+ B \ + % + B i ta k * u k (i), l 

“* s’K + B \ + % + V“k - \ un ) 


(C79) 


where 


B^ = i-th column of B 

u^d) = i-th component of the control vector u^ 

The term B^a^ accounts for the failure and the term B^u k (i) accounts for the loss 
of that control input. 


52 



APPENDIX D 


THE KALMAN FILTER 


As described in appendices A and B, to perform the FDI computations the GLR 
technique uses the innovations from a discrete time Kalman Filter, which estimates 
the system states. As discussed in the Results section of this report, a brief 
investigation was conducted to evaluate the performance of the GLR algorithm as a 
function of the number of states in the filter. From the results of this investiga- 
tion it was decided for this work to have a 10-state Kalman Filter estimate the wind 
gusts as well as the aircraft states. This appendix discusses the implementation of 
this filter. 


The System Model 

The model of the system which the Kalman Filter is designed to estimate must 
include the dynamics of the aircraft and of the wind. This is accomplished by 
combining the aircraft dynamics in equation (C54) and the wind gust dynamics in 
equation (C44) into one state equation as follows: 


V" + * f \ * I 'F S k 

where the combined system state vector is defined by 
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The control vector for the combined system is the same as for the aircraft system and 
is defined as in equation (C3) by 


— . — - 


— —« 

\ (,) 


fie 

\ <2 > 


fith 

\ <3) 


fist 

_\ (4) _ 


fisp 


(D3) 


The plant disturbance vector (winds) is a combination of the white Gaussian noise 
vector sequence and the gust components of the white Gaussian noise sequence 

defined by equations (C56) and (C46), respectively, 


V 1 * 


^( 2 ) 

V 3) 

TfcU) 

V 5 > 

nk (6) 

\d> 

5 k <2) 

5 k <3> 

S k <4) 


(D4) 


The combined system transition matrix is 





0 , <J> 

Y w 


(D5) 


F 

where <J> is the aircraft system transition matrix defined by equation (C39), <}> w is 

the upper left 4x4 sub-matrix of the wind system transition matrix defined by 
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equation (C45), and 1^ is the first four columns of the matrix defined by equa- 
tion (C55). The combined control transition matrix is defined by 




'l* 

0 


D6) 


where i|> is the aircraft control transition matrix defined by equation (C42) . The 
combined disturbance transtion matrix is the identity matrix. 


r 


= I 


(D7 ) 


The plant noise covariance is modeled by 




5 i* 


R 





(D8) 


(D9) 


where R and R' are the covariances of n and respectively. R is 

defined fly equation (C60) , and R' is the upper left 4x4 sub-matrix of the matrix 
defined by equation (C50) . 

For use in the Kalman Filter design the system measurements, or observations, 
are modeled as linear combinations of the system state plus additive Gaussian noise. 


r k+1 


= Hx 


k+1 


+ v: 


k+1 


(DIO) 


The actual measurements, as described in appendix C, are not linear functions of the 
state in all cases. A pre-processor, to be described later, is used to convert the 
measurements z k into a vector y k of observations, or pseudo-measurements, which 
can then be approximated as linear functions of the state. 
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The Filter 


The Kalman Filter is defined by the equations 


x (k+l|k) = ((ipXtklk) + iJ^u^ (Dll) 

x(k|k) = x(k|k-1) + K y(k) (D12) 


where K is the Kalman gain. The innovations y(k) are related to the observations 
by the following: 


y(k) = y k 


Hx(k | k-1 ) 


(D13) 


A constant gain filter was used to greatly reduce the computational workload. 
Constant gain could be used since the aircraft system matrices are constant, or very 
slowly changing, for a given flight segment. The wind covariance may vary more 
rapidly, but since it is also not accurately known, a constant value was used. The 
filter gain K was computed using the subroutine ASYMFIL from the linear control 
system design software package ORACLS (ref. 22). 


The Pre-Processor 

As previously noted, the actual measurements are not all linear functions 

of the state. A pre-processor was used to obtain from the measurements a set of 
observations y^ which could be approximated by a linear model as in equation (DIO). 
This pre-processor is described in the following paragraphs. 

Pitch .- From equation (C62) the output of the pitch attitude gyro is 


z k (1 ) = 0 + v k (1 ) + b(1 ) 


In this case it suffices to let 


Y k (1 ) = « k (1 ) (D14) 

and then y (1) can be modeled as 

y k d) = e + v£(i) 

= x£(1 ) + v£(1 ) (D15) 
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Accelerations The outputs of the body mounted accelerometers are given by 
equations (C63) through (C66). The first step in pre-processing these measurements 
is to convert them to stability coordinates according to 


F F F 

A = z, (2) cos a + z. (3) sin a 
xs k ok o 

(D16) 

F F F 

A = -z, (2) sin a + z, (3) cos a 
zs k ok o 

where a = trim angle-of-attack in the combined system model 

Comparison of the stability axis accelerations in equations (D16) with those in (C65) 
shows that 


P 

A = A + e 
XS XS X 


A 


F 


zs 


A + 
zs 


e 

z 


(D17) 


where e x is an error term due to accelerometer noise Vj c (2), bias error b(2), 
scale factor error s, accelerometer alignment error a, and error (a - a ) in the 
trim angle-of-attack. From equations (C63) and (C64) ° ° 


A xs = U o e 2 ( te k + Bu k + ^k) + qU o tan ° + g sin (0 o + 0 " “o* 
A ZS = “o ^ ° h (^k + BU k + ^k) ‘ qU o" g COS(6 o + 6 “ 


(D18) 




+ u' )U /cos 
o 


«] 


(Ax + Bu. + Dw.) 


Define the observations y k (2) and y k (3) as follows: 


- 7 [<. - « «£- <>l - h b f\ 

o 

y k (3) = ~F «S + g COS (0 o“ a o } ] " e 3 ^ 


(D19) 
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The pre-processor then calculates the observations y k (2) and y k (3) 
measurements z k (2) and z k (3) us i n 9 equations (D16) and (D19). From 
tions (D18) and (D19) these observations can be approximated as linear 
the states as follows: 


y k (2) « q tan a + e^Vk + d f w J 


+ (9/Uq) cos(6^- /) sin 0 + v£(2) 

6 6 

l A p (2,j) x£(j) + l D co (2,j) x^ (j+6) 

j=1 j=1 


+ (g/U F ) cos(0 F - /) X F (1) + V F (2) 

v o ' v o o ' k k 


y k (3) » -q + tan a • e 2 [A^ + D p w k J 


where 


+ [(1 + u’)/cos 2 a] • e^[A ,x + D w ] + vf(3) 

.3 F K. F K K 


F F 

+(g/U ) sinf0 - a ) sin 0 

O v O o ' 


-x k (4) + (g/a') sin(0 F - /) x k (1) 


+ .1 A p ( 3 , j ) x k (j) + l D ( 3 , j ) x F (j+6) + vf(3) 


j=1 


j=1 


D 


co 


F F 
D C 

wo 


from the 
equa- 

functions of 


(D20 ) 


(D21 ) 


(D22) 
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and 


C 

co 



w 


6=0 


-sm a 


-cos a 


/■ F Fx 

-cos (6 - a 1 

V o o ' 


sinfe - a ) 
o o' 


-cos a 


sm a 


• ( F-, 

■sin 6 - a I 

x r> n J 


r F Fx 

-COS (6 - a 1 

^ O O' 


TTV 


TTV 


4b 


4b 


Pitch rate .- From equation (C64) the output of the pitch rate gyro is 
z^(4) = q + v^(4) 

It is sufficient to let 
y k ( . 4) = z k (4) 

which case y^(4)can be modeled as 
y k (4) = x£(4) + \(4) 

Air speed .- From equation (C68) the air speed measurement is 
z k (5) = V a [1 v k (5)] + b(5) 

= U [(1 + u' + u ' ) /cos ( a + a )] • [1 + v, (5)] + b(5) 
o*- w w J k 

Define the air speed observation as 


y k (5) = [z k (5)/U^.] - 1 


(D23) 


(D24) 


(D25) 


(D26) 
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The observation can be approximated as a linear combination of the states as follows: 


y k < 5 ) = 


1 + u' + u' 
w 

cos ( a + a x ) 


[1 + v k (5)] + -^1 - i 


u' + u' + v (5) 
W k 


" x^(2) + l C wQ (1,j) x£(j+6) + v£(5) 

j=1 


Altitude rate.- From equation (C69) the altitude rate measurment is 


z v (6) = h + v v (6) 


= U Q (1 + u' ) [sin( 0 Q + 6 - a 0 ) - tan a cos(0 o + 0 - o^ ) )] + v k (6) 


Define the observation 


y k ( 6 ) - [z k (6)/0^] - (e^ - <£) 


which can be approximately modeled as 


y k (6) “ (1 + u')[sin(0 F - a F ) cos 0 + cos(0* - a*) sin 0 


F F- 


- tan a cos(0 F - a F ] cos 0 + tan a sinfo* - a fc ) sin 0] 

v O O ' O O ' 


F F- 


+ V (6)/U e - (0 F - a F ) 
k o o o' 


u’ sinf0 F - a F ] + cosfo* - 
' o o' ' o 


a F ) 0 - a + v F (6) 
o ' V 


(£- 0\ (2 > ♦<»> -<»> *^(«» 


(D27 ) 


(D28) 


(D29 ) 
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Angle-of-attack . - From equation (C70) the output of the angle-of-attack sensor 
is 


* k (7> - a + + v k (7) + b(7) 


Define the observation 


y k (7) » a k (7) (D30) 


which can be modeled by 


y k <7> - a + a* + v‘<7) 


«x‘(3) * l 0^(2,;}) xJ(j+6) + v£<7) 
j=1 


(D31 ) 


The pre-processor computes the obsrvations y k from the measurements a k 
according to equations (D14), (016), (D19), (D24), (D26), (D28), and (D30). These 
observations . are the inputs to the Kalman Filter. For use in the filter algorithms 
the observations can now be expressed as in equation (DIO) with the observation 
matrix defined by equation (D32). 
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TABLE Ill.b.- GLR VALUES FOR 20-SAMPLE WINDOW 



♦Used 10 ft/sec steady state wind in the x-direction. Other thunderstorm runs used 
39 ft/sec. 
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TABLE III.C.- GLR VALUES FOR 10-SAMPLE WINDOW 


Failure 


Type 


Elevator 


10 c 


Turbulence 

None 

Medium 

Heavy 

Thunder- 

storm 


10180 

9334* 


Throttle 


10870* 

9241* 


Stabilizer 


19320 

17030* 


Spoiler 


30900 

27670* 


17030 


Used 10 ft/sec steady state wind in the x-direction. Other thunderstorm runs used 
39 ft/sec. 






















































TABLE IV.- FULL GLR FAILURE DETECTION PERFORMANCE WITH 
THRESHOLDS SET AT NO TURBULENCE LEVEL 


Failure 

Turbulence 

Type 

Magnitude 

None 

Medium 

Heavy 

Thunder- 

storm 

Elevator 

10° 


Y/Y/Y 

Y/Y/Y 

Y/Y/Y 

Y/Y/Y 

3° 

Y/Y/Y 

Y/Y/Y 

Y/Y/Y 

Y/Y/Y 

Y/Y/Y. 


1° 

Y/Y/n 

n/n/Y 



Throttle 

40° 


. 

Y/Y/Y 

Y/Y/Y 

Y/Y/Y 

Y/Y/Y 

12° 

Y/Y/Y 

Y/Y/Y 

Y/Y/Y 

Y/Y/Y 

Y/Y/Y 


4° 

Y/Y/Y 

Y/Y/Y 



Stabilizer 

-6° 

Y/Y/Y 

Y/Y/Y 

Y/Y/Y 

Y/Y/Y 

Y/Y/Y 

Y/Y/Y 

Y/Y/Y 

3° 

Y/Y/Y 

Y/Y/Y 

Y/Y/Y 


-1° 

Y/Y/Y 

Y/Y/Y 



Spoiler 

8° 

Y/Y/Y 

Y/Y/Y 

Y/Y/Y 

Y/Y/Y 

Y/Y/Y 

3° 

Y/Y/Y 

Y/Y/Y 

Y/Y/Y 

Y/Y/Y 

Y/Y/Y 


1° 

Y/Y/Y 

Y/Y/Y 




Y/Y/Y indicates detection using the 30- , 20- , and 10-sample windows, 
n indicates failure not detected. 























































TABLE V.- FULL GLR FAILURE DETECTION PERFORMANCE WITH 
THRESHOLDS SET AT HEAVY TURBULENCE LEVEL 


Failure 

Turbulence 

Type 

Magnitude 

None 

Medium 

Heavy 

Thunder- 

storm 

Elevator 

10° 


Y/Y/Y 

Y/Y/Y 

Y/Y/Y 

Y/Y/Y 

3° 

n/n/n 

n/n/n 

n/n/n 

n/n/n 

n/n/n 


1° 

n/n/n 

n/n/n 



Throttle 

40° 


Y/Y/Y 

Y/Y/Y 

Y/Y/Y 

Y/Y/Y 

1 29 

Y/n/n 

n/n/n 

Y/n/n 

n/n/n 

Y/Y/Y 


4° 

n/n/n 

n/n/n 



Stabilizer 

-6° 

y/y/y 

Y/Y/Y 

Y/Y/Y 

Y/Y/Y 

Y/Y/Y 

Y/Y/Y 

Y/Y/Y 

3° 

n/n/n 

n/n/n 

n/n/n 


-1° 

n/n/n 

n/n/n 



Spoiler 

8° 

Y/Y/Y 

Y/Y/Y 

Y/Y/Y 

Y/Y/Y 

Y/Y/Y 

3° 

n/n/n 

n/n/n 

n/n/n 

Y/Y/Y 

n/n/n 


1 ° 

n/n/n 

n/n/n 




Y/Y/Y indicates detection using the 30- , 20- , and 10-sample windows, 
n indicates failure not detected. 

































































































TABLE VII.- CONSTRAINED GLR FAILURE DETECTION PERFORMANCE WITH 
THRESHOLDS SET AT HEAVY TURBULENCE LEVEL 


Failure 

Turbulence 

Type 

Magnitude 

None 

Medium 

Heavy 

Thunder- 

storm 

Elevator 

10° 


Y/Y/Y 

Y/Y/Y 

Y/Y/Y 

3° 

Y/Y/Y 

n/Y/Y 

n/Y/Y 

n/Y/Y 
n/Y/Y ■ 


1 ° 

n/n/n 

n/n/n 



Throttle 

40° 


Y/Y/Y 

Y/Y/Y 

Y/Y/Y 

12° 

n/n/n 

n/n/n 

Y/n/n 

n/n/n 

Y/Y/Y 


4° 

n/n/n 

n/n/n 



Stabilizer 

-6° 


Y/Y/Y 

Y/Y/Y 

Y/Y/Y 

3° 

Y/Y/Y 

Y/Y/Y 

Y/Y/Y 

n/n/n 

Y/Y/Y 


-1 ° 

n/n/n 

n/n/n 



Spoiler 

8° 


Y/Y/Y 

Y/Y/Y 

Y/Y/Y 

3° 

Y/n/n 

Y/n/n 

Y/n/n 

Y/Y/Y 

Y/n/n 


1 ° 

n/n/n 

n/n/n 




Y/Y/Y indicates detection using the 30- , 20- , and 10-sample windows, 
n indicates failure not detected. 




TABLE VIII.- CONSTRAINED GLR FAILURE IDENTIFICATION PERFORMANCE: 
FAILURE VECTOR CORRECTLY IDENTIFIED 


Fail 

ure 

Turbulence 

Type 

Magnitude 

None 

Medium 

Heavy 

Thunder- 

storm 

Elevator 

10° 


Y/Y/Y 

Y/Y/Y 

Y/Y/n 

3° 

Y/Y/Y 

*/Y/Y 

*/Y/Y 



1 0 

*/ */ * 

*/ */ * 



Throttle 

O 

o 


Y/Y/Y 

Y/Y/Y 

Y/Y/Y 

12° 

*/ */ * 

*/ */ * 
Y/*/* 

*/ */ * 
Y/Y/Y 


4° 

*/ */ ★ 

*/ */ * 



Stabilizer 

-6° 


Y/Y/n 

Y/Y/n 

n/n/n 

3° 

Y/Y/n 

Y/Y/n 

Y/Y/n 

*/ */ * 
Y/Y/n 


-1 ° 

*/ */ * 

★ / */ * 



Spoiler 

8° 


Y/n/n 

Y/Y/n 

n/n/n 

3° 

Y/*/* 

Y/*/* 

Y/*/* 

n/Y/n 

Y/*/* 


1° 

*/ */ * 

★ / */ * 




Y/Y/Y indicates correct identification using the 30-, 20-, and 10-sample windows, 
n indicates incorrect identification. 

* indicates failure not detected. 
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TABLE IX.- CONSTRAINED GLR FAILURE IDENTIFICATION PERFORMANCE: PERCENT 

ERROR IN FAILURE MAGNITUDE ESTIMATE 


Fai: 

Lure 

Turbulence 

Type 

Magnitude 

None 

Medium 

Heavy 

Thunder- 

storm 

Elevator 

10° 

mmm 

-1 .41# 

-2.21 

-3.40 




-2.58 

-2.30 

2.90 




-0.99 

-0.36 

* 


3° 

-6.93 

* * 

* * 




-9.13 

-9.80 -7.47 

-6.00 -4.43 




-4.20 

-3.77 6.60 

-1.13 14.8 



1 ° 







n 






lll» 

HK9HHI 



Throttle 

40° 

■■ 

44.8 

8.38 

-288.0 




92.0 

54.5 

-485.0 



1 

217.0 

165.0 

-990.0 


12° 

1 

* 88.8 

* 135.0 





* * 

* 214.0 




WEBB 

* * 

* 407.0 


' 

4° 

■ 

■ 











* 

* 




#For each case the percent errors are listed in order of decreasing 
size; that is, the error for the 30-sample window islisted first, 
errors for the 20- and 10-sample windows. 

*Indicates failure not detected. 

^Indicates failure not correctly identified. 


data window 
followed by the 



TABLE IX.- Concluded 


Failure 

Turbulence 

Type 

Magnitude 

None 

Medium 

Heavy 

Thunder- 

storm 

Stabilizer 

-6° 


-37.7# 

-50.2 

k k 




-68.3 

-87.2 

k k 




★ * 

k k 

k k 



4.80 

-3.53 0.97 

* 8.10 




35.4 

26.3 36.0 

* 49.0 




* * 

** ** 

* ** 



-i° 

* 

* 





* 

* 





* 

* 



Spoiler 

8° 


-0.74 

-2.14 

kk 




** 

2.23 

k k 




** 

* ★ 

k k 



22.4 

20.1 27.5 

** 17.6 




* 

•k k 

47.2 * 




* 

k k 

kk k 


« 

i ° 

* 

k 





* 

k 





* 

k 




#For each case the percent errors are listed in order of decreasing data window 
size; that is, the error for the 30-sample window is listed first, followed by the 
errors for the 20- and 10-sample windows. 

^Indicates failure not detected, 
vindicates failure not correctly identified. 
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Figure 3.- Full GLR results with no turbulence and -6° stabili 

failure at 3 seconds. 
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Figure 5.- Histogram of delay time between failure occurrence 
and failure detection for constrained GLR. 
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